% ***************************************************************************************
% m_calibration_spillover.m 
% ***************************************************************************************
clear all; 
close all;
format shortG;
% ==================== SET FUNDAMENTAL PARAMETERS ========================
gravity = 8;             % gravity coefficient (estimation in step2 gravity)   
rho = 0.66;              % discount factor (after 1955) 
rho50 = 0.66;            % --- (for 1950) 
sigma = rho/gravity;     % shape parameter of idiosyncratic shocks (after 1955)  
sigma50 = rho50/gravity; % --- (for 1950)

mu=0.15;                 % expenditure share on floor spaces (consumers) 
gammaL=0.7;              % cost share on labor (production)
gammaH=0.1;              % cost share on floor spaces (production)
eta=4;                   % floor space supply elasticity

% grid of alpha, beta 
lb = -0.25; 
ub = 0.50; 
gridint = 0.001;  
alphagrid=(lb:gridint:ub);  
betagrid=(lb:gridint:ub);      
Ngrid=size(alphagrid,2); 

% **************************************************************
% ******* SET PARAMETERS FOR SPILLOVERS (SPATIAL DECAY) ********
% **************************************************************
spillA_vec = (0.1:0.01:2); 
spillB_vec = (0.1:0.01:2); 
Nj = size(spillA_vec,2); 

%% ========================= LOAD DATA  =========================================  
datafile = '../../data/processed data/quant/Hiroshima_DATA.csv';  CITYDATA = readtable(datafile); 

ID = CITYDATA.geocode;   IDunique = unique(ID); 
if size(ID,1) ~= size(IDunique,1) 
    disp(' >>> CHECK DUPLICATED LOCATION DATA <<<');  pause 
end 
clear IDunique 
N=size(ID,1); Area = CITYDATA.area; GEOGRAPHY=[ID Area]; cbddist=CITYDATA.cbddist; 

Rdata36=CITYDATA.pop36; Rdata45=CITYDATA.pop45bomb; 
Rdata51=CITYDATA.pop51; Rdata55=CITYDATA.pop55; 
Rdata60=CITYDATA.pop60; Rdata65=CITYDATA.pop65;     
Rdata70=CITYDATA.pop70; Rdata75=CITYDATA.pop75; 
Rdata50=CITYDATA.pop50; Rdata53=CITYDATA.pop53;
Rdata57 = 0.6*Rdata55+0.4*Rdata60; 
Rdata66 = 0.8*Rdata65+0.2*Rdata70;

Ldata38=CITYDATA.emp38; Ldata45=CITYDATA.emp45bomb; 
Ldata53=CITYDATA.emp53; Ldata57=CITYDATA.emp57; 
Ldata63=CITYDATA.emp63; Ldata66=CITYDATA.emp66;     
Ldata75=CITYDATA.emp75; 

% Read latitude and longitude (for readability of the csv output in GIS)
lat = CITYDATA.lat; lon = CITYDATA.lon;

Ldata50=Ldata53*(sum(Rdata50)/sum(Rdata53)); 
Ldata55=0.5*Ldata53+0.5*Ldata57;  
Ldata60=0.5.*Ldata57+0.5.*Ldata63; 
Ldata65=(1/3)*Ldata63 + (2/3)*Ldata66; 
Ldata70=(5/9).*Ldata66+(4/9).*Ldata75; 

Rdata45=(sum(Ldata45)/sum(Rdata45)).*Rdata45;  
Rdata50=(sum(Ldata50)/sum(Rdata50)).*Rdata50;
Rdata55=(sum(Ldata55)/sum(Rdata55)).*Rdata55;  
Rdata60=(sum(Ldata60)/sum(Rdata60)).*Rdata60;
Rdata65=(sum(Ldata65)/sum(Rdata65)).*Rdata65;  
Rdata70=(sum(Ldata70)/sum(Rdata70)).*Rdata70;
Rdata75=(sum(Ldata75)/sum(Rdata75)).*Rdata75; 

% data for the post-war economy (every 5 years from 1950);
% NOTE: we use data from 1950 to 1975 for estimation
Rdata=[Rdata50 Rdata55 Rdata60 Rdata65 Rdata70 Rdata75];
Ldata=[Ldata50 Ldata55 Ldata60 Ldata65 Ldata70 Ldata75]; 
totalR = sum(Rdata,1); totalL = sum(Ldata,1);  
Popdata=totalR; 
if min(totalR ~= totalL) == 1   % check total pop ; 
    disp('>>> CHECK TOTAL POPULATION AND EMPLOYMENT ! <<< ');   pause
end 

% all data (every 5 years, 1945-75) 
Rdata_all = [Rdata45 Rdata50 Rdata55 Rdata60 Rdata65 Rdata70 Rdata75];
Ldata_all = [Ldata45 Ldata50 Ldata55 Ldata60 Ldata65 Ldata70 Ldata75]; 

% population/employment density; 
NT=size(Rdata,2);   
Rdensdata = Rdata./repmat(Area,1,NT);  
Ldensdata=Ldata./repmat(Area,1,NT);  

% COMPUTE LOWER BOUNDS for Theta; 
Rthetalb = (Rdata(:,2:NT)-Rdata(:,1:NT-1) < 0);  
Rthetalb = (1 - Rdata(:,2:NT)./Rdata(:,1:NT-1)).*Rthetalb;  
Lthetalb = (Ldata(:,2:NT)-Ldata(:,1:NT-1) < 0);  
Lthetalb = (1 - Ldata(:,2:NT)./Ldata(:,1:NT-1)).*Lthetalb; 
Rthetalb = max(Rthetalb); Lthetalb = max(Lthetalb); 
thetalb = ceil(100.*max([Rthetalb; Lthetalb]))./100;

% ================= LOAD FLOORSPACE PRICES ================================ 
l_q50=CITYDATA.lnQ_50; q50=exp(l_q50*(4/eta)); 
l_q55=CITYDATA.lnQ_55; q55=exp(l_q55*(4/eta)); 
l_q60=CITYDATA.lnQ_60; q60=exp(l_q60*(4/eta)); 
l_q65=CITYDATA.lnQ_65; q65=exp(l_q65*(4/eta)); 
l_q70=CITYDATA.lnQ_70; q70=exp(l_q70*(4/eta));
l_q75=CITYDATA.lnQ_75; q75=exp(l_q75*(4/eta)); 
Qdata=[q50 q55 q60 q65 q70 q75]; % Floor space price matrix  

% ================= LOAD COMMUTING COST ===================================
load('commuting_cost.mat', 'Kappa_mat','TT_walk');  

% ================= SET Theta VECTOR ====================================== 
thetavec=[0.53; 0.53; 0.53; 0.53; 0.53];

if min(thetavec - thetalb')<0
    disp('>>> CHECK THETA VALUES (LOWER BOUND) ! <<< ');   pause
end 

%  =============  GROUP TRACTS FOR GMM BY DISTANCE FROM CBD ===============
%  =============  N=174 BLOCKS to 5 GROUPS (BINS) ==>> 35 x 4 + 34 ======== 
% INDEX LOCATIONS BY ASCENDING ORDER OF CBD distance (FROM SMALL TO LARGE)
% sortn = location index when sorting locations by density from small to large
G=5;
[Gmat] = gmmiv(cbddist,G,N);  GG = sum(Gmat,2); 

%  ========================================================================
%  ==================  MAIN INVERSION OF THE MODEL   ======================
%  ========================================================================
options0 = optimset('Display','none','Algorithm','trust-region-dogleg','MaxFunEvals',50000000, 'MaxIter',1000000,'TolFun',1e-6,'TolX',1e-6);
Xmat=zeros(N,NT); 
Ymat=zeros(N,NT); 
Kmat = zeros(N,N,NT); 
LQmat=zeros(N,NT); 

%  construct kernel for previous periods NT-1, ..., 1;  
Kmat(:,:,NT)=Kappa_mat(:,:,NT).^(-rho/sigma);   % last period 
for t=1:NT-1
    tt=NT-t; 
    theta=thetavec(tt); 
    if t < NT-1; Kmat(:,:,tt)=(Kappa_mat(:,:,tt).^(-rho/sigma)).*(Kmat(:,:,tt+1).^(rho*(1-theta))); 
    else;        Kmat(:,:,tt)=(Kappa_mat(:,:,tt).^(-rho50/sigma50)).*(Kmat(:,:,tt+1).^(sigma*rho50*(1-theta)/sigma50)); end
end 

% CONSTRUCT CONTINUATION VALUES FOR FLOOR SPACE PRICES 
LQmat(:,NT)=Qdata(:,NT);
for t=1:NT-1 
    tt=NT-t;
    theta=thetavec(tt); 
    if t < NT-1; LQmat(:,tt)=Qdata(:,tt).*(LQmat(:,tt+1).^(rho*(1-theta))); 
    else;        LQmat(:,tt)=Qdata(:,tt).*(LQmat(:,tt+1).^(rho50*(1-theta))); end
end 

K50=Kmat(:,:,1);
Q50=LQmat(:,1); 

%  =================== INVERSION OF FUNDAMENTALS =========================== 
for t=1:NT-1
    tt = NT-t;    % compute from the last period (1975,1970,1965,1960,1955)
    theta=thetavec(tt); 
    K = Kmat(:,:,tt+1); 
    Q=LQmat(:,tt+1);
    Rpre=Rdata(:,tt); Rpost=Rdata(:,tt+1); 
    Lpre=Ldata(:,tt); Lpost=Ldata(:,tt+1);
       
    X0=ones(N,1);  Y0=ones(N,1);     
    Rsyst = @(x)fReval2(x,K,Q,Rpre,Rpost,Lpre,Lpost,rho,sigma,theta,mu,gammaL,gammaH);
    Lsyst = @(y)fLeval2(y,K,Q,Rpre,Rpost,Lpre,Lpost,rho,sigma,theta,mu,gammaL,gammaH);  
    [x_fsolve, xfval_fsolve]=fsolve(Rsyst, X0, options0); X=abs(x_fsolve); X=X./geomean(X);  
    [y_fsolve, yfval_fsolve]=fsolve(Lsyst, Y0, options0); Y=abs(y_fsolve); Y=Y./geomean(Y);
    Xmat(:,tt+1)=X;    Ymat(:,tt+1)=Y; 
end 
X55 = Xmat(:,2); Y55 = Ymat(:,2); 
X60 = Xmat(:,3); Y60 = Ymat(:,3); 
X65 = Xmat(:,4); Y65 = Ymat(:,4); 
X70 = Xmat(:,5); Y70 = Ymat(:,5); 
X75 = Xmat(:,6); Y75 = Ymat(:,6);

estimateA_mat = zeros(Nj,2); 
estimateB_mat = zeros(Nj,2); 

%  ===================== DECOMPOSE FUNDAMENTALS ===========================
for jj = 1:Nj

spillA = spillA_vec(jj); 
spillB = spillB_vec(jj);
spillA_mat = exp(-spillA.*TT_walk); 
spillB_mat = exp(-spillB.*TT_walk); 

amat3D=zeros(N,NT-1,Ngrid); 
bmat3D=zeros(N,NT-1,Ngrid);
a_errmat3D = zeros(N,NT-1,Ngrid); 
b_errmat3D=zeros(N,NT-1,Ngrid);
for ii=1:Ngrid
  
    alpha = alphagrid(ii); beta = betagrid(ii);
    [a,a_err] = ffund_spill(alpha,Ymat,Ldensdata,rho,sigma,spillA_mat,thetavec,N,NT);     
    [b,b_err] = ffund_spill(beta,Xmat,Rdensdata,rho,sigma,spillB_mat,thetavec,N,NT);

    amat3D(:,:,ii)=a;  a_errmat3D(:,:,ii)=a_err; 
    bmat3D(:,:,ii)=b;  b_errmat3D(:,:,ii)=b_err; 

end 

%  ===========  GRID SEARCH: 1st Step GMM OBJECTIVE FUNCTION ===============
objvala = zeros(Ngrid,1);  
objvalb = zeros(Ngrid,1);
d_a_err = a_errmat3D(:,2:NT-1,:)-a_errmat3D(:,1:NT-2,:);   
d_b_err = b_errmat3D(:,2:NT-1,:)-b_errmat3D(:,1:NT-2,:);

for ii=1:Ngrid
    objvala(ii)=fobjm(Gmat,GG,d_a_err(:,:,ii),N,NT,G,eye(G*(NT-2))); 
    objvalb(ii)=fobjm(Gmat,GG,d_b_err(:,:,ii),N,NT,G,eye(G*(NT-2)));
end

[~, ii_mina]  = min(objvala);
[~, ii_minb]  = min(objvalb);

alphaest1st = alphagrid(ii_mina);
betaest1st = betagrid(ii_minb); 
[alphaest1st, betaest1st]

% ======================================================================
% ==========================  2nd Step GMM   ===========================
% ======================================================================
d_a_err1st = d_a_err(:,:,ii_mina);
d_b_err1st = d_b_err(:,:,ii_minb); 

mmadist1st=zeros(G,N,NT-2); mmbdist1st=zeros(G,N,NT-2);
for t=1:NT-2
    mmadist1st(:,:,t)=Gmat.*repmat(d_a_err1st(:,t)',G,1);  mmbdist1st(:,:,t)=Gmat.*repmat(d_b_err1st(:,t)',G,1);   
end
va = zeros(G*(NT-2),G*(NT-2)); vb=va; 
for i=1:N
    mmadistv = reshape(squeeze(mmadist1st(:,i,:)), G*(NT-2),1);  
    mmbdistv = reshape(squeeze(mmbdist1st(:,i,:)), G*(NT-2),1);  
    va = va + mmadistv*mmadistv'./N;   vb = vb + mmbdistv*mmbdistv'./N; 
end

objvala2=zeros(Ngrid,1);  
objvalb2=zeros(Ngrid,1);
wa = inv(va); wb=inv(vb);    % weight matrix 
for ii=1:Ngrid
    objvala2(ii)=fobjm(Gmat,GG,d_a_err(:,:,ii),N,NT,G,wa);
    objvalb2(ii)=fobjm(Gmat,GG,d_b_err(:,:,ii),N,NT,G,wb); 
end

[~, ii_mina2]  = min(objvala2);
[~, ii_minb2]  = min(objvalb2);

alphaest = alphagrid(ii_mina2);
betaest = betagrid(ii_minb2); 
[alphaest betaest] 

estimateA_mat(jj,1) = spillA; 
estimateA_mat(jj,2) = alphaest; 

estimateB_mat(jj,1) = spillB; 
estimateB_mat(jj,2) = betaest; 

end 


%% image of objective function 
f=gcf;
set(f,'PaperPositionMode','auto'); set(f,'PaperOrientation','landscape');   
set(f,'Position', [0, 0, 800, 800]);   
xx = estimateA_mat(:,1); yy = estimateA_mat(:,2); 
xx2 = estimateB_mat(:,1); yy2 = estimateB_mat(:,2); 
plot(xx,yy,'k','LineWidth',3); hold on
plot(xx2,yy2,'--b','LineWidth',3); hold off
set(gca,'FontSize',13)
xlim([0.2 2.0])
ylim([0.05 0.30])
legend('Productivity spillovers', 'Amenity spillovers','Fontsize', 14);
xlabel('Spatial decay in productivity and amenity spillovers','interpreter','latex','fontsize',17,'fontweight','bold');
ylabel('Agglomeration parameter ($\alpha$, $\beta$)','interpreter','latex','fontsize',17,'fontweight','bold');
print(gcf, '-dpdf', '../../output/figure/figure_spillover_agglomeration.pdf');
close all; 

%% ============================= END ======================================= 

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%                       FUNCTIONS
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ---------- FUNCTION FOR CONSTRUCTING IV --------- 
function[gmat]=gmmiv(ivvar,nG,N); 
% ivvar: variable for constructing IV
Ground=round(N/nG); gn=zeros(N,1); gmat=zeros(nG,N);
[~,sortn] = sort(ivvar,'ascend'); 
Gthres=[Ground+1; 2*Ground+1; 3*Ground+1; 4*Ground+1; 5*Ground+1; 6*Ground+1; 7*Ground+1; 8*Ground+1; 9*Ground+1];
for i=1:N
    ii=find(sortn==i);
    if ii<Gthres(1); gn(i)=1; 
    elseif ii >= Gthres(nG-1); gn(i)=nG; 
    else; 
        for g=2:nG-1 
        if Gthres(g-1) <= ii && ii < Gthres(g); gn(i)=g;  
        else; continue  
        end 
        end 
    end 
end
for i=1:N
    gi = gn(i);
    gmat(gi,i)=1;
end
end 

% --------  FUNCTION FOR COMPUTING CHARACTERISTICS BASED ON GRAVITY STURCTURE --------
function[Reval] = fReval2(X,K,Q,DRpre,DRpost,DLpre,DLpost,rho,sigma,theta,mu,gammaL,gammaH)   % eq D.16. X corresponds to Theta

N = size(DRpre,1); 
X=abs(X);  
X=X./geomean(X);
xtemp=K.*((repmat(Q',N,1).^(-mu)).*repmat(X',N,1)).^(rho/sigma);
xtemp=xtemp./repmat(sum(xtemp,2),1,N); 

Rdif=DRpost-(1-theta).*DRpre; 
Ldif=DLpost-(1-theta).*DLpre; 

Reval = Rdif - sum(xtemp.*repmat(Ldif,1,N),1)'; 
Reval = abs(Reval); 

end

function[Leval] = fLeval2(Y,K,Q,DRpre,DRpost,DLpre,DLpost,rho,sigma,theta,mu,gammaL,gammaH) % eq D.17. Y corersponds to Omega

N = size(DRpre,1); 
Y=abs(Y);  
Y=Y./geomean(Y);
ytemp= K.*((repmat(Q,1,N).^(-gammaH/gammaL)).*(repmat(Y,1,N).^(1/gammaL))).^(rho/sigma);
ytemp=ytemp./repmat(sum(ytemp,1),N,1); 

Rdif=DRpost-(1-theta).*DRpre; 
Ldif=DLpost-(1-theta).*DLpre; 

Leval = Ldif - sum(ytemp.*repmat(Rdif',N,1),2); 
Leval = abs(Leval); 

end

% --------  FUNCTION FOR COMPUTING THE EXOGENOUS FUNDAMENTALS --------
function[c, c_err] = ffund_spill(param,Zmat,densdata,rho,sigma,spill_mat,thetavec,N,NT)  

Z0=Zmat(:,NT); cmat=zeros(N,NT-1);
for t=1:NT-1
    denspost=densdata(:,NT-t+1); 
    spillover=spill_mat*denspost; 
    invspill = spillover.^(-param);
    cc = invspill.*Z0;
    cmat(:,NT-t)= cc; 
    theta=thetavec(NT-t); 
    Z0=Zmat(:,NT-t).*(Z0.^(-rho*(1-theta)));     
end

c = log(cmat);   % N x NT-1
c_te = sum(c,1)./N;
c_err = c-repmat(c_te,N,1); 
end 

% ----------  FUNCTION FOR OBJECTIVE FUNCTION ------------
% Gmat: G (number of grids) x N (number of locations)  
% cmat: N (number of locations) x NT-2 (number of periods well defined for changes)
% Wmat: weight matrix, G*NT-2 (number of grids) x G*NT-2 
% objdist: scalar 
function[objdist] = fobjm(Gmat,GG,cmat,N,NT,G,Wmat); 
mdist = zeros(G,N,NT-2);
for t=1:NT-2
    mmt = Gmat.*repmat(cmat(:,t)',G,1);  % G x N
    mdist(:,:,t) = mmt;                 
end
mdist = squeeze(sum(mdist,2));   %  G x NT-2
mdist = reshape(mdist, G*(NT-2),1); 
objdist=(mdist./N)'*Wmat*(mdist./N);
end
